! written by jxzou at 20220220: generate .47 file from Gaussian .fch(k) file

program main
 use util_wrapper, only: formchk
 implicit none
 integer :: i
 character(len=240) :: fchname

 i = iargc()
 if(i /= 1) then
  write(6,'(/,A)') ' ERROR in program fch247: wrong command line argument!'
  write(6,'(A,/)') ' Example: fch247 a.fch'
  stop
 end if

 call getarg(1, fchname)
 call require_file_exist(fchname)

 ! if .chk file provided, convert into .fch file automatically
 i = LEN_TRIM(fchname)
 if(fchname(i-3:i) == '.chk') then
  call formchk(fchname)
  fchname = fchname(1:i-3)//'fch'
 end if

 call fch247(fchname)
end program main

subroutine fch247(fchname)
 use fch_content
 implicit none
 integer :: i, j, k, m, nhigh, fid
 integer, allocatable :: eff_nuc_charge(:), cen_assign(:)
 real(kind=8), allocatable :: S(:,:),F(:,:),dm_a(:,:),dm_b(:,:), D(:,:,:)
 character(len=1), parameter :: l_str(0:5) = ['S','P','D','F','G','H']
 character(len=240) :: f47
 character(len=240), intent(in) :: fchname
 logical :: uhf

 call check_uhf_in_fch(fchname, uhf)
 call read_fch(fchname, uhf)

 i = INDEX(fchname, '.fch')
 f47 = fchname(1:i-1)//'.47'
 open(newunit=fid,file=TRIM(f47),status='replace')
 write(fid,'(A,I0,A,I0)',advance='no') ' $GENNBO  NATOMS=',natom,'  NBAS=',nbf
 if(uhf) write(fid,'(A)',advance='no') '  OPEN'
 write(fid,'(A)') '  UPPER  BODM  $END'
 write(fid,'(A)') ' $NBO  $END'
 write(fid,'(A)') ' $COORD'
 write(fid,'(A)') ' Generated by utility fch247 in MOKIT'

 ! calculate effective nuclear chages (remember to minus ECP core, if any)
 allocate(eff_nuc_charge(natom), source=0)
 if(allocated(RNFroz)) then
  forall(i = 1:natom) eff_nuc_charge(i) = ielem(i) - INT(RNFroz(i))
 else
  forall(i = 1:natom) eff_nuc_charge(i) = ielem(i)
 end if
 do i = 1, natom, 1
  write(fid,'(1X,2I6,3F15.6)') ielem(i),eff_nuc_charge(i),coor(:,i)
 end do ! for i
 deallocate(eff_nuc_charge)
 write(fid,'(A)') ' $END'
 write(fid,'(A)') ' $BASIS'

 allocate(cen_assign(nbf))
 k = 0
 do i = 1, ncontr, 1
  j = shltyp2nbf(shell_type(i))
  cen_assign(k+1:k+j) = shell2atom_map(i)
  k = k + j
 end do ! for i
 deallocate(shell2atom_map)
 write(fid,'(A,17I7)') '  CENTER =',cen_assign(1:min(17,nbf))
 if(nbf > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,nbf)

 ! Next, we use the array cen_assign to store LABEL
 cen_assign = 0; k = 0
 do i = 1, ncontr, 1
  j = shell_type(i)
  select case(j)
  case(0) ! S
   cen_assign(k+1) = 1
  case(1) ! P
   forall(m = 1:3) cen_assign(k+m) = 100 + m
  case(2) ! 6D
   cen_assign(k+1:k+6) = [201,204,206,202,203,205]
  case(3) ! 10F
   cen_assign(k+1:k+10) = [301,307,310,304,302,303,306,309,308,305]
  case(4) ! 15G
   forall(m = 1:15) cen_assign(k+m) = 416 - m
  case(5) ! 21H
   forall(m = 1:21) cen_assign(k+m) = 522 - m
  case(-1) ! L/SP
   cen_assign(k+1:k+4) = [1,101,102,103]
  case(-2) ! 5D
   cen_assign(k+1:k+5) = [255,252,253,254,251]
  case(-3) ! 7F
   forall(m = 1:7) cen_assign(k+m) = 350 + m
  case(-4) ! 9G
   forall(m = 1:9) cen_assign(k+m) = 450 + m
  case(-5) ! 11H
   forall(m = 1:11) cen_assign(k+m) = 550 + m
  case default
   stop
  end select
  k = k + shltyp2nbf(j)
 end do ! for i

 write(fid,'(A,17I7)') '   LABEL =',cen_assign(1:min(17,nbf))
 if(nbf > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,nbf)
 deallocate(cen_assign)
 write(fid,'(A)') ' $END'
 write(fid,'(A)') ' $CONTRACT'
 write(fid,'(A,I7)') '  NSHELL =',ncontr
 write(fid,'(A,I7)') '    NEXP =',nprim

 allocate(cen_assign(ncontr))
 forall(i = 1:ncontr) cen_assign(i) = shltyp2nbf(shell_type(i))
 write(fid,'(A,17I7)') '   NCOMP =',cen_assign(1:min(17,ncontr))
 if(ncontr > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,ncontr)

 write(fid,'(A,17I7)') '   NPRIM =',prim_per_shell(1:min(17,ncontr))
 if(ncontr > 17) write(fid,'(10X,17I7)') (prim_per_shell(i),i=18,ncontr)

 cen_assign(1) = 1
 do i = 2, ncontr, 1
  cen_assign(i) = cen_assign(i-1) + prim_per_shell(i-1)
 end do ! for i

 write(fid,'(A,17I7)') '    NPTR =',cen_assign(1:min(17,ncontr))
 if(ncontr > 17) write(fid,'(10X,17I7)') (cen_assign(i),i=18,ncontr)
 deallocate(cen_assign)

 write(fid,'(A,4E16.7)') '     EXP =',prim_exp(1:min(4,nprim))
 if(nprim > 4) write(fid,'(10X,4E16.7)') (prim_exp(i),i=5,nprim)

 ! (1) multiply contraction coeff by normalization factor
 ! (2) use F(:,:) to store CS, CP, CD, CF, CG, CH.
 i = MAXVAL(shell_type)
 if(i < 0) i = -i
 j = MINVAL(shell_type)
 if(j < 0) j = -j
 nhigh = MAX(i,j)
 allocate(F(nprim,0:nhigh),source=0d0)
 m = 0
 do i = 1, ncontr, 1
  j = prim_per_shell(i); k = shell_type(i)
  if(k == -1) then ! L, SP
   call contr_coeff_multiply_norm_fac(j,0,prim_exp(m+1:m+j),contr_coeff(m+1:m+j))
   call contr_coeff_multiply_norm_fac(j,1,prim_exp(m+1:m+j),contr_coeff_sp(m+1:m+j))
   F(m+1:m+j,0) = contr_coeff(m+1:m+j)
   F(m+1:m+j,1) = contr_coeff_sp(m+1:m+j)
  else !  S/P/D/F/G/H
   call contr_coeff_multiply_norm_fac(j,k,prim_exp(m+1:m+j),contr_coeff(m+1:m+j))
   F(m+1:m+j,ABS(k)) = contr_coeff(m+1:m+j)
  end if
  m = m + j
 end do ! for i
 do i = 0, nhigh, 1
  write(fid,'(A,4E16.7)') '      C'//l_str(i)//' =',F(1:min(4,nprim),i)
  if(nprim > 4) write(fid,'(10X,4E16.7)') (F(j,i),j=5,nprim)
 end do ! for i
 deallocate(F)
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $OVERLAP'
 allocate(S(nbf,nbf))
 if(nbf == nif) then
  call solve_ovlp_from_cct(nbf, alpha_coeff, S)
 else
  write(6,'(/,A)') 'ERROR in subroutine fch247: nbf/=nif, linear dependence&
                  & detected.'
  write(6,'(A)') 'Solving overlap integrals not supported currently.'
  stop
 end if
 write(fid,'(2X,5E15.7)') ((S(j,i),j=1,i),i=1,nbf)
 deallocate(S)
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $DENSITY'
 if(uhf) then ! unretricted type wave function
  if(allocated(tot_dm) .and. allocated(spin_dm)) then
   dm_a = 0.5d0*(tot_dm + spin_dm) ! alpha density
   dm_b = 0.5d0*(tot_dm - spin_dm) ! beta density
  else
   allocate(dm_a(nbf,nbf), dm_b(nbf,nbf))
   dm_a = MATMUL(alpha_coeff(:,1:na), TRANSPOSE(alpha_coeff(:,1:na)))
   dm_b = MATMUL(beta_coeff(:,1:nb), TRANSPOSE(beta_coeff(:,1:nb)))
  end if
  write(fid,'(2X,5E15.7)') ((dm_a(j,i),j=1,i),i=1,nbf)
  write(fid,'(2X,5E15.7)') ((dm_b(j,i),j=1,i),i=1,nbf)
  deallocate(dm_a, dm_b)
 else ! retricted type wave function
  if(.not. allocated(tot_dm)) then
   allocate(tot_dm(nbf,nbf))
   tot_dm = 2d0*MATMUL(alpha_coeff(:,1:nb), TRANSPOSE(alpha_coeff(:,1:nb)))
   if(na > nb) tot_dm = tot_dm + MATMUL(alpha_coeff(:,nb+1:na), &
                                        TRANSPOSE(alpha_coeff(:,1:nb+1:na)))
  end if
  write(fid,'(2X,5E15.7)') ((tot_dm(j,i),j=1,i),i=1,nbf)
 end if
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $FOCK'
 allocate(F(nbf,nbf))
 call solve_fock_from_ctfc(nbf, nif, alpha_coeff, eigen_e_a, F)
 write(fid,'(2X,5E15.7)') ((F(j,i),j=1,i),i=1,nbf)
 if(uhf) then
  call solve_fock_from_ctfc(nbf, nif, beta_coeff, eigen_e_b, F)
  write(fid,'(2X,5E15.7)') ((F(j,i),j=1,i),i=1,nbf)
 end if
 deallocate(F)
 write(fid,'(A)') ' $END'

 write(fid,'(A)') ' $LCAOMO'
 write(fid,'(2X,5E15.7)') alpha_coeff
 if(uhf) write(fid,'(2X,5E15.7)') beta_coeff
 write(fid,'(A)') ' $END'

! write(fid,'(A)') ' $DIPOLE'
! allocate(D(nbf,nbf,3))
! call calc_dipole_mat_using_fch(fchname, nbf, D)
! do i = 1, 3
!  write(fid,'(2X,5E15.7)') ((D(k,j,i),k=1,j),j=1,nbf)
! end do ! for i
! deallocate(D)
! write(fid,'(A)') ' $END'

 close(fid)
end subroutine fch247

! multiply each contraction coefficient with corresponding normalization factor
! Note: this subroutine is different with that in file fch2wfn.f90
subroutine contr_coeff_multiply_norm_fac(n, p, prim_exp, contr_coeff)
 implicit none
 integer :: i, p0
 integer, intent(in) :: n, p
 real(kind=8) :: r2a_pi(n)
 ! 2*alpha/PI. Because variable cannot begin with a integer, so 'r' is added
 real(kind=8), parameter :: PI = 4d0*DATAN(1d0)
 real(kind=8), parameter :: sqrt_2pi = DSQRT(8d0*DATAN(1d0))
 real(kind=8), intent(in) :: prim_exp(n)
 real(kind=8), intent(inout) :: contr_coeff(n)

 if(p == -1) then
  write(6,'(A)') 'ERROR in subroutine contr_coeff_multiply_norm_fac: p=-1.'
  write(6,'(A)') 'You should divide L/SP into separate S/P before calling &
                    &this subroutine.'
  stop
 end if
 if(p < 0) then
  p0 = -p
 else
  p0 = p
 end if

 forall(i = 1:n)
  r2a_pi(i) = (2d0*prim_exp(i)/PI)**(0.25d0)
  contr_coeff(i) = contr_coeff(i)*(r2a_pi(i)**(2*p0+3))
 end forall
 if(p0 > 0) contr_coeff = contr_coeff*(sqrt_2pi**p0)

 select case(p0)
 case(0,1) ! do nothing
 case(2) ! D
  contr_coeff = contr_coeff/DSQRT(3d0)
 case(3) ! F
  contr_coeff = contr_coeff/DSQRT(15d0)
 case(4) ! G
  contr_coeff = contr_coeff/DSQRT(105d0)
 case(5) ! H
  contr_coeff = contr_coeff/DSQRT(945d0)
 case default
  write(6,'(A,I0)') 'ERROR in subroutine contr_coeff_multiply_norm_fac: inva&
                       &lid p0=', p0
  stop
 end select
end subroutine contr_coeff_multiply_norm_fac

